#include "StdAfx.h"
#include "Ambigous_Zone_Based_Extraction.h"
#include "LinkedList.h"


Ambigous_Zone_Based_Extraction::Ambigous_Zone_Based_Extraction(void)
{
	this->source_image = NULL;
	this->half_thinned_image = NULL;
	this->thinned_image = NULL;
	this->Sf.clear();
	this->Se.clear();
	this->S1.clear();
	this->Sa.clear();
}

int Ambigous_Zone_Based_Extraction::run(const char* filename, vector<StrokeNode*> &stroke_graph, IplImage** ambigous_zones_map)
{
	this->source_image = cvLoadImage(filename,CV_LOAD_IMAGE_GRAYSCALE);
	if (!this->source_image)
	{
		cerr << "Couldn't open the image file" << endl;
		return 1;
	}
	
	this->half_thinned_image = cvCreateImage(cvSize(this->source_image->width, this->source_image->height), this->source_image->depth, this->source_image->nChannels);
	this->thinned_image = cvCreateImage(cvSize(this->source_image->width, this->source_image->height), this->source_image->depth, this->source_image->nChannels);
	cvCopy(this->source_image, this->half_thinned_image);
	cvCopy(this->source_image, this->thinned_image);

	Half_thinning_Algorithm *thinning_algorithm = new Half_thinning_Algorithm();
	this->w = thinning_algorithm->calculate_estimated_stroke_width(source_image);
	thinning_algorithm->run_half_thinning(this->half_thinned_image);
	cvSaveImage("half_thinned_image.png",this->half_thinned_image);
	thinning_algorithm->run(this->thinned_image);
	delete thinning_algorithm;

	this->find_Candidate_Fork_Points();
	this->find_Spurious_Candidate_Fork_Points();
	this->visualize_fork_points();
	this->process_Fork_Points();
	
	IplImage* img = this->detect_Ambigous_Zone(stroke_graph);
	*ambigous_zones_map = cvCloneImage(img);
	cvReleaseImage(&img);
	return 0;
}

void Ambigous_Zone_Based_Extraction::find_Candidate_Fork_Points()
{
	// Create mask
	int** mask = new int * [3];
			
	for (int k = -1; k < 2; k++)
	{
		mask[k+1] = new int [3];
	}

	for (int y = 0; y< (int) this->thinned_image->height; y++)
	{
		for (int x = 0; x< (int) this->thinned_image->width; x++)
		{
			if (((uchar *)(this->thinned_image->imageData + y * this->thinned_image->widthStep))[x] == 0)
			{

				// Fill mask values
				for (int k = -1; k < 2; k++)
				{
					for (int l = -1; l < 2; l++)
					{
						if ((y + k < 0) || (x + l < 0) || (y + k >= this->thinned_image->height) || (x + l >= this->thinned_image->width) )
						{
							mask[k+1][l+1] = 0;
						}
						else
						{
							uchar pixel_value = (((uchar *)(this->thinned_image->imageData + (y + k) * this->thinned_image->widthStep))[x + l] == 255) ? 0 : 1;
							mask[k+1][l+1] = pixel_value;
						}
					}
				}

				int Nc = 0;
				int Nb = 0;
				int ordered_set[9] = {mask[0][1],mask[0][2],mask[1][2],mask[2][2],mask[2][1],mask[2][0],mask[1][0],mask[0][0], mask[0][1]};
				
				for (int i = 0; i< 8; i++)
				{
					if ((ordered_set[i] == 0) && (ordered_set[i+1] == 1))
					{
						Nc++;
					}
					
					Nb+= ordered_set[i];
				}
				
				if (Nb == 1)
				{
					CvPoint p;
					p.x = x;
					p.y = y;
					this->Se.insert(p);
				}
				else if ((Nc >= 3) || (Nb >= 4))
				{
					CvPoint p;
					p.x = x;
					p.y = y;
					this->Sf.insert(p);
				}


			}
		}
	}

	// Delete mask
	for (int k = 0; k < 3; k++)
	{
		delete [] mask[k];
	}
	delete [] mask;	

	return;
}

 

void Ambigous_Zone_Based_Extraction::find_Spurious_Candidate_Fork_Points()
{
	set<CvPoint,ImagePointComp>::iterator it;
	set<CvPoint,ImagePointComp> branch;
	for (it = this->Se.begin(); it != this->Se.end(); it++)
	{
		branch.clear();
		vector<CvPoint> childs;
		childs.clear();
		childs.push_back(*it);
		while(!childs.empty())
		{
			CvPoint currentNode = childs.back();
			childs.pop_back();
			branch.insert(currentNode);
			if (this->Sf.find(currentNode) != this->Sf.end())
			{
				double Lb = sqrt((((*it).x - currentNode.x) * ((*it).x - currentNode.x) + ((*it).y - currentNode.y) * ((*it).y - currentNode.y)) * 1.0);
				double d1 = 0.8 * this->w;
				double dh = 2.5 * this->w;
				if (Lb < d1)
				{
					this->S1.insert(currentNode);
					for (int k = -1; k < 2; k++)
					{
						for (int l = -1; l < 2; l++)
						{
							if ((currentNode.y + k < 0) || (currentNode.x + l < 0) || (currentNode.y + k >= this->source_image->height) || (currentNode.x + l >= this->source_image->width) )
							{
								continue;
							}
							else
							{
								if (((uchar *)(this->thinned_image->imageData + (currentNode.y + l) * this->thinned_image->widthStep))[currentNode.x + k] == 0)
								{
									CvPoint temp;
									temp.x = currentNode.x + k;
									temp.y = currentNode.y + l;
									if (branch.find(temp) == branch.end())
									{
										if (this->Sf.find(temp) != this->Sf.end())
										{
											this->S1.insert(temp);
										}
									}
								}
							}
						}
					}

				}
				else
				{
					if (Lb <= dh)
					{
						double sum = 0.0;
						int Np = (int) branch.size();
						double max_dc = 0.0;
						set<CvPoint,ImagePointComp>::iterator it_temp;
						for (it_temp = branch.begin(); it_temp != branch.end(); it_temp++)
						{
							double dc = this->calculate_dc((*it_temp));
							sum+= dc;
							if (max_dc < dc)
							{
								max_dc = dc;
							}
						}
						sum = sum / Np;
						if ((sum > 0.6 * this->w) || (max_dc > this->w))
						{
							this->S1.insert(currentNode);
							for (int k = -1; k < 2; k++)
							{
								for (int l = -1; l < 2; l++)
								{
									if ((currentNode.y + k < 0) || (currentNode.x + l < 0) || (currentNode.y + k >= this->source_image->height) || (currentNode.x + l >= this->source_image->width) )
									{
										continue;
									}
									else
									{
										if (((uchar *)(this->thinned_image->imageData + (currentNode.y + l) * this->thinned_image->widthStep))[currentNode.x + k] == 0)
										{
											CvPoint temp;
											temp.x = currentNode.x + k;
											temp.y = currentNode.y + l;
											if (branch.find(temp) == branch.end())
											{
												if (this->Sf.find(temp) != this->Sf.end())
												{
													this->S1.insert(temp);
												}
											}
										}
									}
								}
							}
						}
					}
				}
				//break;
			}
			else
			{
				for (int k = -1; k < 2; k++)
				{
					for (int l = -1; l < 2; l++)
					{
						if ((currentNode.y + k < 0) || (currentNode.x + l < 0) || (currentNode.y + k >= this->source_image->height) || (currentNode.x + l >= this->source_image->width) )
						{
							continue;
						}
						else
						{
							if (((uchar *)(this->thinned_image->imageData + (currentNode.y + l) * this->thinned_image->widthStep))[currentNode.x + k] == 0)
							{
								CvPoint temp;
								temp.x = currentNode.x + k;
								temp.y = currentNode.y + l;
								if (branch.find(temp) == branch.end())
								{
									childs.push_back(temp);
								}
							}
						}
					}
				}
			}
		}
	}
	
	return;
}

void Ambigous_Zone_Based_Extraction::process_Fork_Points()
{
	set< CvPoint,ImagePointComp > S2;

	// S2 = Sf - S1;
	set< CvPoint,ImagePointComp >::iterator it;
	for (it = this->Sf.begin(); it != this->Sf.end(); it++)
	{

		if (this->S1.find(*it) == this->S1.end())
		{
			S2.insert(*it);
		}
	}

	// Create S4 set
	set< CvPoint,ImagePointComp > S4;
	for (int y = 0; y< (int) this->thinned_image->height; y++)
	{
		for (int x = 0; x< (int) this->thinned_image->width; x++)
		{
			if (((uchar *)(this->thinned_image->imageData + y * this->thinned_image->widthStep))[x] == 0)
			{
				CvPoint pi;
				pi.x = x;
				pi.y = y;
				S4.insert(pi);
			}
		}
	}

	// Create S3 set
	set<CvPoint, ImagePointComp> S3;
	S3 = Tools::marking_borders(this->half_thinned_image);

	// Create (S4 - S3) set
	set<CvPoint, ImagePointComp> S43;
	for (it = S4.begin(); it != S4.end(); it++)
	{
		if (S3.find(*it) == S3.end())
		{
			S43.insert(*it);
		}
	}

	// Create S5 = S2 ^ (S4 - S3)
	set<CvPoint, ImagePointComp> S5;
	for (it = S2.begin(); it != S2.end(); it++)
	{
		if (S43.find(*it) != S43.end())
		{
			S5.insert(*it);
		}
	}

	Graph *fork_point_graph = new Graph(S2, this->thinned_image);
	set< CvPoint,ImagePointComp > S2_clone = S2;

	while (!S43.empty())
	{
		vector<CvPoint> childs;
		childs.push_back(*S43.begin());
		S43.erase(*S43.begin());
		set<CvPoint,ImagePointComp> e;
		while (!childs.empty())
		{
			CvPoint currentPoint = childs.back();
			childs.pop_back();
			if (S2_clone.find(currentPoint) != S2_clone.end())
			{
				e.insert(currentPoint);
				S2_clone.erase(currentPoint);
			}

			for (int y = -1; y < 2; y++)
			{
				for (int x = -1; x < 2; x++)
				{
					if ((x != 0) || (y != 0))
					{
						CvPoint neighbourPoint;
						neighbourPoint.x = currentPoint.x + x;
						neighbourPoint.y = currentPoint.y + y;
						if (S43.find(neighbourPoint) != S43.end())
						{
							childs.push_back(neighbourPoint);
							S43.erase(neighbourPoint);
						}
					}
				}
			}
		}
		if (e.size() > 1)
		{
			this->Sa.push_back(e);
			for (set< CvPoint,ImagePointComp >::iterator i = e.begin(); i!= e.end(); i++)
			{
				S2.erase((*i));
			}
		}
	}

	for (int i = 0; i< (int) this->Sa.size(); i++)
	{
		set<CvPoint,ImagePointComp>::iterator it;
		for (it = this->Sa[i].begin(); it != this->Sa[i].end(); it++)
		{
			CvPoint currentPoint = (*it);
			vector<CvPoint> branchPoints = fork_point_graph->get_branch_points(currentPoint.x, currentPoint.y);
			for (int j = 0; j< (int) branchPoints.size(); j++)
			{
				if (S2.find(branchPoints[j]) != S2.end())
				{
					double distance = sqrt((double)((currentPoint.x - branchPoints[j].x) * (currentPoint.x - branchPoints[j].x) + (currentPoint.y - branchPoints[j].y) * (currentPoint.y - branchPoints[j].y)));

					if (distance < 1.4 * this->w)
					{
						this->Sa[i].insert(branchPoints[j]);
						S2.erase(branchPoints[j]);
					}
				}
			}
		}
	}



	//while (!S5.empty())
	//{
	//	ImagePoint currentPoint = *S5.begin();
	//	set<ImagePoint,ImagePointComp> e;
	//	e.insert(currentPoint);
	//	
	//	vector<ImagePoint> branchPoints = fork_point_graph->get_branch_points(currentPoint.x, currentPoint.y);
	//	for (int j = 0; j< (int) branchPoints.size(); j++)
	//	{
	//		if ((S5.find(branchPoints[j]) == S5.end()) && (e.find(branchPoints[j]) == e.end()))
	//		{

	//			double distance = sqrt((double)((currentPoint.x - branchPoints[j].x) * (currentPoint.x - branchPoints[j].x) + (currentPoint.y - branchPoints[j].y) * (currentPoint.y - branchPoints[j].y)));
	//			if (distance < 0.8 * this->w)
	//			{
	//				e.insert(branchPoints[j]);
	//				S2.erase(branchPoints[j]);
	//				cout << "(" << branchPoints[j].x << "," << branchPoints[j].y << "), ";
	//			}
	//		}
	//	}

	//	S5.erase(currentPoint);
	//	S2.erase(currentPoint);
	//	cout << "{(" << currentPoint.x << "," << currentPoint.y << "), ";
	//	vector<ImagePoint> neighbours = fork_point_graph->find_neighbours(currentPoint.x, currentPoint.y);

	//	for (int i = 0; i< (int) neighbours.size(); i++)
	//	{
	//		if (S5.find(neighbours[i]) != S5.end())
	//		{
	//			e.insert(neighbours[i]);
	//			vector<ImagePoint> branchPoints = fork_point_graph->get_branch_points(neighbours[i].x, neighbours[i].y);
	//			for (int j = 0; j< (int) branchPoints.size(); j++)
	//			{
	//				if ((S5.find(branchPoints[j]) == S5.end()) && (e.find(branchPoints[j]) == e.end()))
	//				{

	//					double distance = sqrt((double)((neighbours[i].x - branchPoints[j].x) * (neighbours[i].x - branchPoints[j].x) + (neighbours[i].y - branchPoints[j].y) * (neighbours[i].y - branchPoints[j].y)));
	//					if (distance < 0.8 * this->w)
	//					{
	//						e.insert(branchPoints[j]);
	//						S2.erase(branchPoints[j]);
	//						cout << "(" << branchPoints[j].x << "," << branchPoints[j].y << "), ";
	//					}
	//				}
	//			}
	//			S5.erase(neighbours[i]);
	//			S2.erase(neighbours[i]);
	//			cout << "(" << neighbours[i].x << "," << neighbours[i].y << "), ";
	//		}
	//	}
	//	cout << "}" << endl;
	//	this->Sa.push_back(e);
	//}

	while (!S2.empty())
	{
		CvPoint currentPoint = *S2.begin();
		set<CvPoint,ImagePointComp> e;
		e.insert(currentPoint);
		S2.erase(currentPoint);

		vector<CvPoint> branchPoints = fork_point_graph->get_branch_points(currentPoint.x, currentPoint.y);
		for (int j = 0; j< (int) branchPoints.size(); j++)
		{
			if (S2.find(branchPoints[j]) != S2.end())
			{

				double distance = sqrt((double)((currentPoint.x - branchPoints[j].x) * (currentPoint.x - branchPoints[j].x) + (currentPoint.y - branchPoints[j].y) * (currentPoint.y - branchPoints[j].y)));

				if (distance < 1.4 * this->w)
				{
					e.insert(branchPoints[j]);
					S2.erase(branchPoints[j]);
				}
			}
		}
		this->Sa.push_back(e);
		
	}

	//for (int i = 0; i< (int) this->Sa.size(); i++)
	//{
	//	set<ImagePoint,ImagePointComp>::iterator it;
	//	cout << "{";
	//	for (it = this->Sa[i].begin(); it != this->Sa[i].end(); it++)
	//	{
	//		ImagePoint currentPoint = (*it);
	//		cout << "(" << currentPoint.x << "," << currentPoint.y << "), ";
	//	}
	//	cout << "}" << endl;
	//}

	delete fork_point_graph;
}

double Ambigous_Zone_Based_Extraction::calculate_dc(CvPoint pi)
{
	double dc;
	Half_thinning_Algorithm *thinning_algorithm = new Half_thinning_Algorithm();
	dc = thinning_algorithm->calculate_dc(pi.x, pi.y, this->source_image);
	delete thinning_algorithm;
	return dc;
}
void Ambigous_Zone_Based_Extraction::visualize_fork_points()
{
	IplImage *display_image = cvCreateImage(cvSize(this->source_image->width, this->source_image->height), IPL_DEPTH_8U, 3);
	for (int y = 0; y< (int) this->thinned_image->height; y++)
	{
		for (int x = 0; x< (int) this->thinned_image->width; x++)
		{
			uchar pixel_value = ((uchar *)(this->thinned_image->imageData + y * this->thinned_image->widthStep))[x];
			((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 0] = pixel_value;
			((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 1] = pixel_value;
			((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 2] = pixel_value;
		}
	}

	
	set<CvPoint,ImagePointComp>::iterator it;
	for (it = this->Sf.begin(); it != this->Sf.end(); it++)
	{
		int x = (*it).x;
		int y = (*it).y;

		((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 0] = 0;
		((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 1] = 0;
		((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 2] = 255;
	}

	for (it = this->Se.begin(); it != this->Se.end(); it++)
	{
		int x = (*it).x;
		int y = (*it).y;
		((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 0] = 255;
		((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 1] = 0;
		((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 2] = 0;
	}

	for (it = this->S1.begin(); it != this->S1.end(); it++)
	{
		int x = (*it).x;
		int y = (*it).y;
		((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 0] = 44;
		((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 1] = 189;
		((uchar *)(display_image->imageData + y * display_image->widthStep))[x * 3 + 2] = 28;
	}
	cvSaveImage("fork_points.png",display_image);
	cvReleaseImage(&display_image);
	return;
}

IplImage* Ambigous_Zone_Based_Extraction::detect_Ambigous_Zone(vector<StrokeNode*> &stroke_graph)
{
	stroke_graph.clear();

	IplImage *img_edge = cvCreateImage(cvSize(this->source_image->width,this->source_image->height),this->source_image->depth,this->source_image->nChannels);


	CvMemStorage* storage = cvCreateMemStorage(); 
	
	CvChain* first_contour = NULL;

	// Invert the image 
	cvThreshold(this->source_image, img_edge, 100, 255, CV_THRESH_BINARY_INV);

	int Nc = cvFindContours(
		img_edge,
		storage,
		(CvSeq**)&first_contour,
		sizeof(CvChain),
		CV_RETR_LIST,
		CV_CHAIN_CODE);
	
	vector<LinkedList*> contour_list;
	contour_list.clear();

	for( CvChain* c=first_contour; c!=NULL; c=(CvChain*)c->h_next )
	{
		LinkedList *contour = new LinkedList();

		CvChainPtReader reader;
		int i, total = c->total;
		if (c->total > 0)
		{
			cvStartReadChainPoints( c, &reader);
			for( i = 0; i < total; i++ )
			{
				CvPoint point;
				CV_READ_CHAIN_POINT(point, reader);
				int code = reader.code;
				contour->insert(point, code);
			}
			contour_list.push_back(contour);
		}
		
	}

	cvReleaseImage(&img_edge);
	for( CvChain* c=first_contour; c!=NULL; c=(CvChain*)c->h_next )
	{
		cvClearSeq((CvSeq*)c);
	}
	cvClearMemStorage(storage);
	cvReleaseMemStorage(&storage);

	IplImage* des = cvCreateImage(cvSize(this->source_image->width, this->source_image->height),8,3);
	cvCvtColor(this->source_image, des,CV_GRAY2BGR);
	

	// For each element e in Sa
	for (int i = 0; i< (int) this->Sa.size(); i++)
	{
		// Build the Sc set
		set<CvPoint,ImagePointComp>::iterator it;
		CvPoint centerPoint;
		centerPoint.x = 0;
		centerPoint.y = 0;
		int colorindex = i;

		for (it = this->Sa[i].begin(); it != this->Sa[i].end(); it++)
		{
			CvPoint currentPoint = (*it);

			centerPoint.x+= currentPoint.x;
			centerPoint.y+= currentPoint.y;
		}
		centerPoint.x = centerPoint.x / (int) this->Sa[i].size();
		centerPoint.y = centerPoint.y / (int) this->Sa[i].size();

		// Build Sc set:
		// Sc = {c | d(pc, c) < dr}
		// dr = 1.4 w
		double dr =  2.0 * this->w;
		
		vector<NodeList*> Sc;
		// Iterate all contour points and retrieve relevant points
		for (int i = 0; i< (int) contour_list.size(); i++)
		{
			LinkedList* currentContour = contour_list[i];
			// Traverse all points in current contour
			NodeList* currentPoint = currentContour->get(0);
			
			double d = Tools::distance_Euclidean(centerPoint, currentPoint->data);
			if (d < dr)
			{
				Sc.push_back(currentPoint);
			}

			while(currentPoint->index < currentContour->getLength() - 1)
			{
				currentPoint = currentPoint->next;
				double d = Tools::distance_Euclidean(centerPoint, currentPoint->data);
				if (d < dr)
				{
					Sc.push_back(currentPoint);
				}
			}

		}
		
		// Sz is the set that contains polygons that defines ambigous zone
		vector<NodeList*> Sz;
		
		while (!Sc.empty())
		{

			// Find pm
			// dm = 1.35w
			double dm = (double) 1.2 * this->w;
			double dmin = 10000.0;
			NodeList* pm = NULL;

			for (int i = 0; i< (int) Sc.size(); i++)
			{
				NodeList* currentPoint = Sc[i];
				double d = Tools::distance_Euclidean(centerPoint, currentPoint->data);
				if (dmin > d)
				{
					dmin = d;
					pm = currentPoint;
				}
			}

			if (dmin < dm)
			{
				Sz.push_back(pm);
				int Lth = 0;
				double sin_theta = 0.0;
				double cos_theta = 0.0;
				if (dmin <= this->w * 1.0 / 2)
				{
					sin_theta = 1.0;
					cos_theta = 0.0;
				}
				else
				{
					sin_theta = (this->w * 1.0) / (2 * dmin);
					cos_theta = -1 * sqrt(4.0 * dmin * dmin - this->w * this->w * 1.0) / (2 * dmin);
				}

				Lth = (int) floor(dmin * cos_theta + sqrt(dm * dm - dmin * dmin * sin_theta * sin_theta) + 0.5) ;

				vector<NodeList*> CS;
				if (2 * Lth + 1 >= pm->owner->getLength())
				{
					pm->owner->getSegment(0,pm->owner->getLength() - 1,CS);
				}
				else
				{
					int startIndex = pm->index - Lth;
					if (startIndex < 0)
					{
						startIndex = startIndex + pm->owner->getLength();
					}
					int endIndex = (pm->index + Lth) % pm->owner->getLength();
					pm->owner->getSegment(startIndex,endIndex,CS);
				}

				for (int i = 0; i< (int) CS.size(); i++)
				{
					vector<NodeList*>::iterator it;
					it = find(Sc.begin(), Sc.end(), CS[i]);
					if ( it != Sc.end())
					{
						Sc.erase(it);
					}
				}


			}
			else
			{
				break;
			}
		}


		double Cth = 0.05;
		bool corner_condition = false;
		if (Sz.size() == 3)
		{
			NodeList* v0 = Sz[0];
			NodeList* v1 = Sz[1];
			NodeList* v2 = Sz[2];

			vector<NodeList*> CS;
			double LCS = 0.0;
			
			CS.clear();
			if (v0->owner == v1->owner)
			{
				if (!corner_condition)
				{
					// Retrieve contour segment between v0 and v1
					v0->owner->getSegment(v0->index, v1->index,CS);
					LCS = CS.size();
					// Calculate d0
					double d0 = 0.0;
					for (int i = 0; i< (int) CS.size(); i++)
					{
						double distance = Tools::distance_point_to_line(CS[i]->data,v0->data,v1->data);
						if (d0 < distance)
						{
							d0 = distance;
						}
					}

					double d1 = Tools::distance_Euclidean(v0->data, v1->data) / 2;

					double Larc = 2*(d0 + (d1* d1) / d0) * atan(d0 / d1);

					corner_condition = (abs(LCS - Larc) / Larc) < Cth;

				}
			}
			
			CS.clear();
			if (v0->owner == v2->owner)
			{
				if (!corner_condition)
				{

					// Retrieve contour segment between v0 and v2
					v0->owner->getSegment(v0->index, v2->index,CS);
					LCS = CS.size();
					// Calculate d0
					double d0 = 0.0;
					for (int i = 0; i< (int) CS.size(); i++)
					{
						double distance = Tools::distance_point_to_line(CS[i]->data,v0->data,v2->data);
						if (d0 < distance)
						{
							d0 = distance;
						}
					}

					double d1 = Tools::distance_Euclidean(v0->data, v2->data) / 2;

					double Larc = 2*(d0 + (d1* d1) / d0) * atan(d0 / d1);

					corner_condition = (abs(LCS - Larc) / Larc) < Cth;

				}
			}
			
			CS.clear();
			if (v1->owner == v2->owner)
			{
				if (!corner_condition)
				{

					// Retrieve contour segment between v0 and v2
					v1->owner->getSegment(v1->index, v2->index,CS);
					LCS = CS.size();
					// Calculate d0
					double d0 = 0.0;
					for (int i = 0; i< (int) CS.size(); i++)
					{
						double distance = Tools::distance_point_to_line(CS[i]->data,v1->data,v2->data);
						if (d0 < distance)
						{
							d0 = distance;
						}
					}

					double d1 = Tools::distance_Euclidean(v1->data, v2->data) / 2;

					double Larc = 2*(d0 + (d1* d1) / d0) * atan(d0 / d1);

					corner_condition = (abs(LCS - Larc) / Larc) < Cth;

				}
			}

		}

		
		if ((Sz.size() >= 3) && !corner_condition)
		{
			unsigned int hehe = (unsigned int) ((colorindex + 1) * 50);
			unsigned char r;
			unsigned char g;
			unsigned char b;
			r = (unsigned char) (((unsigned int)(hehe << 24)) >> 24);
			g = (unsigned char) (((unsigned int)(hehe << 16)) >> 24);
			b = (unsigned char) (((unsigned int)(hehe << 8)) >> 24);
			CvScalar green = cvScalar(b * 1.0,g * 1.0,r* 1.0);
			
			vector<CvPoint> points;
			points.clear();
			for (int i = 0; i< (int) Sz.size(); i++)
			{
				points.push_back(Sz[i]->data);
			}
			
			CvPoint* curves[1];
			CvPoint* result;
			vector<CvPoint> zone_polygon_points;
			int nCurvePts[1];
			GrahamImplementation graham_algo;
			result = graham_algo.compute(points, nCurvePts[0], zone_polygon_points);
			curves[0] = result;
			cvFillPoly(des,curves,nCurvePts,1,green);

			StrokeNode* ambigous_zone_node = new StrokeNode();
			
			ambigous_zone_node->elements = Tools::get_substroke_elements(result[0],des,green);
			ambigous_zone_node->polygon_vertex = zone_polygon_points;
			ambigous_zone_node->isAmbigousZone = true;
			ambigous_zone_node->id = colorindex;
			stroke_graph.push_back(ambigous_zone_node);
			delete [] result;
		}
	}	 

	for (int i = 0; i< (int) contour_list.size(); i++)
	{
		delete contour_list[i];
	} 
	return des;
}

Ambigous_Zone_Based_Extraction::~Ambigous_Zone_Based_Extraction(void)
{
	cvReleaseImage(&this->source_image);
	cvReleaseImage(&this->half_thinned_image);
	cvReleaseImage(&this->thinned_image);
}
